library(tidyverse)
library(zoo)
library(cowplot)

# Read ACLED 

acled_mex <- read_csv("ACLED Data_2026-02-06.csv")

# Read VPI

load("fig-5.RData")



acled_mex_by_day <- 
  acled_mex %>%
  dplyr::group_by(event_date) %>%
  dplyr::summarise(tot_fatalities = sum(fatalities)) %>%
  dplyr::mutate(
    fatalities_ma30 = zoo::rollmean(tot_fatalities, k = 30, fill = NA, align = "right")
  )


# Step 1: Calculate 2023 baseline (same as before)
baseline_2023 <- 
  acled_mex_by_day %>%
  filter(event_date >= as.Date("2023-01-01") & event_date <= as.Date("2023-12-31")) %>%
  summarise(acled_baseline = mean(tot_fatalities, na.rm = TRUE)) %>%
  bind_cols(
    daily_range %>%
      filter(date >= as.Date("2023-01-01") & date <= as.Date("2023-12-31")) %>%
      summarise(vpi_baseline = mean(mean_score, na.rm = TRUE))
  )

# Step 2: Full period data (2020-01-01 to 2024-05-31)
temporal_full <- 
  acled_mex_by_day %>%
  rename(date = event_date) %>%
  full_join(
    daily_range %>% select(date, vpi_score = mean_score),
    by = "date"
  ) %>%
  filter(date >= as.Date("2020-01-01") & date <= as.Date("2024-05-31")) %>%
  arrange(date) %>%
  mutate(
    # 30-day moving averages
    acled_ma30 = zoo::rollmean(tot_fatalities, k = 30, fill = NA, align = "center"),
    vpi_ma30 = zoo::rollmean(vpi_score, k = 30, fill = NA, align = "center", na.pad = TRUE),
    
    # Express as fraction of 2023 baseline
    acled_index = acled_ma30 / baseline_2023$acled_baseline,
    vpi_index = vpi_ma30 / baseline_2023$vpi_baseline
  )

# Step 3: Zoom period data
temporal_zoom_1 <- temporal_full %>%
  filter(date >= as.Date("2020-04-15") & date <= as.Date("2020-08-31"))

temporal_zoom_2 <- temporal_full %>%
  filter(date >= as.Date("2023-12-15") & date <= as.Date("2024-05-31"))

# Top panel - Full period
p_top <- temporal_full %>%
  ggplot(aes(x = date)) +
  geom_hline(yintercept = 1, linetype = "dashed", color = "gray50", linewidth = 0.8) +
  geom_line(aes(y = acled_index, color = "ACLED fatalities"), linewidth = 1) +
  geom_line(aes(y = vpi_index, color = "VPI"), linewidth = 1) +
  scale_color_manual(
    values = c("ACLED fatalities" = "#1f77b4", "VPI" = "#d62728"),
    name = ""
  ) +
  scale_x_date(date_breaks = "6 months", date_labels = "%b\n%Y") +
  scale_y_continuous(
    labels = scales::percent_format(accuracy = 1),
    breaks = seq(0, 2, 0.25)
  ) +
  labs(
    y = "Fraction of 2023 average",
    x = NULL,
    title = "Violence indicators relative to 2023 baseline",
    subtitle = "30-day moving averages (2023 average = 100%)"
  ) +
  theme_minimal() +
  theme(
    legend.position = "top",
    panel.grid.minor = element_blank(),
    text = element_text(size = 11),
    plot.title = element_text(size = 13, face = "bold")
  )

# Bottom panel - Zoom period
p_mid <- temporal_zoom_1 %>%
  ggplot(aes(x = date)) +
  geom_hline(yintercept = 1, linetype = "dashed", color = "gray50", linewidth = 0.8) +
  annotate("text", x = as.Date("2020-05-05"), y = 1.05, 
           label = "2023 average", color = "gray50", size = 3.5) +
  geom_line(aes(y = acled_index, color = "ACLED fatalities"), linewidth = 1.3) +
  geom_line(aes(y = vpi_index, color = "VPI"), linewidth = 1.3) +
  scale_color_manual(
    values = c("ACLED fatalities" = "#1f77b4", "VPI" = "#d62728"),
    name = ""
  ) +
  scale_x_date(date_breaks = "1 month", date_labels = "%b\n%Y") +
  scale_y_continuous(
    labels = scales::percent_format(accuracy = 1),
    breaks = seq(0, 2, 0.25)
  ) +
  labs(
    y = "Fraction of 2023 average",
    x = NULL,
    subtitle = "Zoom: April 2020 - August 2020"
  ) +
  theme_minimal() +
  theme(
    legend.position = "none",  # Legend only on top panel
    panel.grid.minor = element_blank(),
    text = element_text(size = 11),
    plot.subtitle = element_text(size = 11, face = "italic")
  )

# Bottom panel - Zoom period
p_bottom <- temporal_zoom_2 %>%
  ggplot(aes(x = date)) +
  geom_hline(yintercept = 1, linetype = "dashed", color = "gray50", linewidth = 0.8) +
  annotate("text", x = as.Date("2024-01-05"), y = 1.05, 
           label = "2023 average", color = "gray50", size = 3.5) +
  geom_line(aes(y = acled_index, color = "ACLED fatalities"), linewidth = 1.3) +
  geom_line(aes(y = vpi_index, color = "VPI"), linewidth = 1.3) +
  scale_color_manual(
    values = c("ACLED fatalities" = "#1f77b4", "VPI" = "#d62728"),
    name = ""
  ) +
  scale_x_date(date_breaks = "1 month", date_labels = "%b\n%Y") +
  scale_y_continuous(
    labels = scales::percent_format(accuracy = 1),
    breaks = seq(0, 2, 0.25)
  ) +
  labs(
    y = "Fraction of 2023 average",
    x = NULL,
    subtitle = "Zoom: December 2023 - May 2024"
  ) +
  theme_minimal() +
  theme(
    legend.position = "none",  # Legend only on top panel
    panel.grid.minor = element_blank(),
    text = element_text(size = 11),
    plot.subtitle = element_text(size = 11, face = "italic")
  )

combined_plot <- plot_grid(
  p_top, 
  p_mid, 
  p_bottom,
  ncol = 1,
  rel_heights = c(1, 1, 1),
  align = "v",
  axis = "lr"  # Align left and right axes
)

# Save if needed
ggsave("vpi_acled_two_panel.jpg", combined_plot, 
       width = 10, height = 15, dpi = 300)


